Non-affine heterogeneities and droplet fluctuations in an equilibrium crystalline solid 
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We show, using molecular dynamics simulations, that a two-dimensional Lennard-Jones solid 
exhibits droplet fluctuations characterized by non-affine deviations from local crystallinity. The 
fraction of particles in these droplets increases as the mean density of the solid decreases and 
approaches « 20% of the total number in the vicinity of the fluid-solid phase boundary. We monitor 
the geometry, local equation of state, density correlations and van Hove functions of these droplets. 
We provide evidence that these non-affine heterogeneities should be interpreted as being droplet 
fluctuations from nearby, metastable minima. The local excess pressure of the droplets plotted 
against the local number density shows a van der Waal's loop with distinct branches corresponding 
to fluid-like compact, and string-like glassy droplets. The distinction between fluid-like and glassy 
droplets disappears above a well defined temperature. 



I. INTRODUCTION 

Crystalline solids typically exhibit local non-affine de- 
formations when driven by external stresses. In many 
instances these non-affine deformations can be described 
in terms of a density of dislocations, however such a de- 
scription is problematic when the density of dislocations 
is large enough that their cores overlap [I]. Since the 
overlapping cores of the dislocations have the character 
of a fluid, it has been suggested that these excitations 
should be thought of as fluid-like droplets [2HS]. This 
has proved a useful interpretation, especially since amor- 
phous solids, for which dislocations are difficult to define, 
also show such localized deformations under shear. 

Just as dislocations in a solid can be thermally excited 
in the absence of external drive, it is reasonable to ask 
whether these fluid-like droplet fluctuations [S] can arise 
in the absence of external perturbation, especially when 
close to the fluid-solid phase boundary. Droplet fluctu- 
ations have been studied in great detail for simple Ising 
systems undergoing a first order transition where they are 
known to influence the asymptotic behavior of dynamic 
correlations and introduce subtle essential singularities 
in the equilibrium free energy [7] ■ The nature and role of 
droplet fluctuations in solids, on the other hand, has not 
received similar attention. In this paper, using a molecu- 
lar dynamics (MD) simulation of a two-dimensional (2D) 
Lennard-Jones (LJ) solid, we show that, indeed, ther- 
mally excited droplet fluctuations do exist close to melt- 
ing. We characterize the local droplet fluctuations using 
a non-affine order parameter [10] . and further classify 
them as being fluid-like or "glassy" (reflecting the whole 
family of non-crystalline metastable configurations). 

We report MD simulations of a 2D single-component 
system with the atoms interacting via LJ potential, viz. 



To = (ma 2 /e) 1 / 2 sets the scale for time with m as the 
mass of the particles. We use e = a = t = 1 without 
loss of generality. The phase diagram of this system as 
obtained from an earlier Monte Carlo study jS] is shown 
in Fig[l] in the scaled temperature T - number density p 
plane. First order liquid-solid and gas-solid boundaries 
are shown. In FigfTJ we have also shown the p and T 
values at which we have obtained our results from equi- 
librated configurations. All our state points lie in the 
single phase region where one always obtains an equi- 
librium, triangular solid. Data showing any evidence of 
local melting of the solid is discarded. 
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where e and a set the scale of energy and length whereas 



FIG. 1: Phase diagram of 2D Lennard-Jones solid as given 
in [8]. The first order boundaries are shown by solid lines. 
Open circles indicate the T and p values at which we have 
performed MD simulations. 



The main results of this paper are summarized below: 

1. We show that there is a significant fraction of non- 
affine droplets in a 2D solid as the density is re- 
duced, the fraction of particles in droplets reaches 
to about 20% at melting. 

2. The droplets are characterized by a density and 
excess pressure over the solid, with positive excess 



pressures associated with string-like droplets and 
negative, with compact droplets. 

3. The density fluctuations of the droplets obeys a dis- 
tinct fluctuation-response relation associated with 
the susceptibility of the droplet. 

4. The excess pressure of a droplet of a given size, as 
we move across the phase diagram, depends on the 
shape of the droplets and is a non-monotonic func- 
tion of the density at low temperatures. At high 
temperatures, this increases monotonically with 
density. 

5. Finally, we show that the equal and unequal time 
density correlations within the droplets are liquid- 
like for the compact and glassy for the string-like 
droplets. 

Taken together, these results suggest that the non-affine 
droplets should be viewed as fluctuations arising from 
nearby metastable liquid and glassy minima. 

The rest of the paper is organized as follows. In the 
next section we give details of our MD simulations and 
the data analysis scheme we use to identify non-afnne 
droplets. We next describe our results for the droplet 
shape, local thermodynamics and density correlations. 
Finally, we discuss the significance of our results and con- 
clude. 



II. SIMULATION AND DATA ANALYSIS 

Our MD simulations are carried out both in the canon- 
ical NVT and micro-canonical NVE ensembles using a 
velocity Verlet algorithm § with a time-step of 10~ 4 . 
Starting from a system of 10 4 particles arranged in a reg- 
ular triangular lattice at desired p — N/V, we have cho- 
sen the initial velocity of each particle from a Maxwell- 
Boltzmann distribution at temperature, T. We equili- 
brate the system at T for an initial 2 x 10 5 MD time 
steps. We then switch to a constant NVE ensemble and 
collect data for another 10 5 MD time steps, storing con- 
figurations at regular intervals. For our system, fluctua- 
tions of T are of order 1 in 10 -4 . At a fixed T and p we 
analyze configurations of particles using a local measure 
for non-affineness (x) defined as the residual deformation 
of a region Q surrounding a particle that is left over after 
fitting the best affine strain measured with respect to the 
ideal triangular lattice at T, p[10]. The neighborhood fi, 
defined using a cutoff distance A and consisting of n par- 
ticles, centered around any tagged particle at r in the 
initial configuration is compared with that of the same 
particle at time t. We obtain the local strain e^ which 
maps as nearly as possible all the n particles from the 
initial to the instantaneous configuration at t. This is 




FIG. 2: (a) Probability distribution of the non-afnne param- 
eter x f° r T = 0.4 at several densities. For these values of 
p, the P(x) are bimodal, the contribution to the first peak 
is shown on the left and for the second peak on the right 
(note change of scale) (b) Comparison of P(x) with that of 
a harmonic solid Pharmfx)- (Inset) Plot of Pd rop = P(x) f° r 
X > Xo and otherwise. 



done by minimizing the (positive) scalar quantity, 
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with respect to e^. Here the indices iandj = x, y and 
r % n {t) and r l n (0) are the i th component of the position 
vector of the n th particle in the initial and instanta- 
neous configurations respectively. Any residual value of 
Xn(r,t), which has units of c 2 , is a measure of non- 
affineness. We have chosen A = 2.5 as our coarse grain- 
ing length. We compute the probability distribution 
P(x) (Fig|2] (a) and (b)) of the coarse grained x- 

Explicit calculations show that x 1S large near de- 
fects such as vacancies and dislocations which result in 
a change in local coordination. For p > 1.0, where the 
solid is expected to be almost harmonic, P(x) shows a 
single peak for \ < 1.0[11|- As density of the system is 
decreased this peak becomes shorter and broader and a 
second peak emerges for 2.5 < x < 5.0. This second peak 
becomes more prominent as the system approaches the 
liquid-solid phase boundary. Inspection of the configu- 
rations of the particles with large x values contributing 
to the second peak in P(x) show that these typically 
represent changes in the local topology where a pair of 
particles from the next near neighbor shell becomes closer 
than their nearest neighbors thereby increasing the local 
density. 



To identify the truly anharmonic droplet fluctuations 
at given p and T, we need to subtract out contribu- 
tions to x coming from purely harmonic distortions. 
In order to do this we note that for a harmonic solid, 
P{x) is unimodal and has a scaling form P(%; T, p) — 
P(xk/A 2 ksTp) where k is the spring constant of the 
harmonic solid. The spring constant of the reference har- 
monic solid was chosen to match the probability distri- 
bution of the lattice parameter in the LJ solid, obtained 
from the curvature of the first peak in the (angle aver- 
aged) pair distribution function g(r) of the solid. The 
probability distribution of x °f the equivalent harmonic 
solid Pharm(x)i was multiplied by a constant till the area 
of the curve matched the area of the first peak of P(x). 
This procedure yields a threshold xo above which there 
is no non-affineness in the reference harmonic solid and 
therefore any non-affiness in the LJ solid above this value 
must necessarily be attributed to anharmonic fluctua- 
tions. This subtraction scheme results in the distribution 
-Pdrop(x) of purely equilibrium droplet fluctuations in a 
solid at constant p and T (Fig Mb) (inset)). 



III. RESULTS 

We have used the above definition of non-affine par- 
ticles to carry out the rest of our analysis. The frac- 
tion of particles in droplets N c /N is typically small but 
increases with decreasing average density, reaching ap- 
proximately 20% close to melting[2 (see Figpl). This 
might seem too large at first, however note that this is 
consistent with typical dislocation densities in 2d solids 
close to melting [15,(12]. We expect this droplet fraction 
to be much lower in 3 dimensions. 
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FIG. 3: Plot of the fraction N c /N of particles in non-affine 
droplets as a function of p. The open symbols corresponds 
to all the non-affine particles based on the threshold criterion 
described in Fig. [3Fb) while the filled symbols correspond 
only to particles with value of \ within the second peak in 
-Pdrop(x)- 



A. Droplet size and shape 

The droplets have a distribution of sizes, shapes, den- 
sity and internal pressure; we compute these quanti- 
ties using standard cluster counting techniques and lo- 
cal Delaunay analysis |12j. The number of particles in 
the droplets n c is exponentially distributed with a mean 
which increases towards the liquid-solid phase boundary. 

In FigEFa) we have shown a typical snapshot of the 
LJ particles for p — 0.92 and T = 0.4. To eliminate 
unimportant random fluctuations we show only droplets 
with n c > 7. The snapshots show both compact and 
string-like morphologies. The droplets are dynamic, they 
coalesce and dissociate while continuously fluctuating in 
shape and size. 
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FIG. 4: (a) Snapshot of a typical particle configuration in 
the LJ solid at T = 0.4 and p = 0.92. The black, filled, cir- 
cles denote non-affine particles while the rest of the particles 
are dots. A string-like and a compact droplet are pointed 
out using red and blue arrows respectively. Close up of a 
(b) string-like and (c) compact droplet showing triangulated 
neighborhoods. While defining droplets, we also include the 
nearest neighbor shell around every non-affine particle. This 
spatial coarse-graining improves statistics. 



To understand shape and size fluctuations of the 
droplets we need to obtain local densities and pressure. 
We do this by constructing a local Delaunay net[H] of 
nearest neighbor particles(see FigEJb) and (c)). The area 
of the droplet A c is then the sum of the areas of the De- 
launay triangles and the density p c = n c /A c . To obtain 
the internal pressure of the droplet, p c , we compute the 
virial (Fy • ry) where Fy and ry are the nearest neigh- 
bor forces and distances respectively for particles i and j 
belonging to the droplet. The average (...) is over all the 
particles n c in the droplet. The droplets are characterized 
by a distribution of p c and excess pressures, Ap c = p c ~p, 
where p is the mean pressure of the surrounding solid. 
Figure [5]ja) shows a scatter diagram of the excess pres- 
sure versus the density in the cluster. We find that the 
droplets with high p c (hence large x) and large (positive) 
Ap c are string-like, whereas droplets with low p c and 
Ap c < are compact (FigMa)). To study the behaviour 



between these two extremes, we argue that these drops 
resemble 2D lattice animals with excess pressure [13j [14] . 
This analogy suggests that the mean radius of gyration 
R g (n c ,Ap c ,T) obeys a crossover relation, 
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where p = Ap c a 2 /ksT. Note that our sign convention 
for the pressure implies that p > corresponds to de- 
flated droplets. Further note that in [13] the number of 
boundary particles are kept fixed which is a different en- 
sembles from ours where the total number of particles in 
the droplets n c is fixed. 

The crossover scaling function asymptotes to, 
F(x— >0) = const, and F(x -> ±oo) = x e± . Note that 
the scaling form takes into account the natural scaling 
p ~ A~ 1 ~ R~ 2 where A c is the area of the droplet. 
At p = we expect that the boundary of the droplet is 
a self avoiding random walk and hence v — 3/4. The 
exponents 

(string-like) and R 



density and pressure of the droplets. In a bulk solid, lo- 
cal thermodynamic equilibrium demands that the local 
variations in the density are related to the pressure com- 
puted from the variation dF/dp of the Hclmholtz free 
energy F with respect to the density via the equation of 
state (EOS) of the solid at the ambient temperature. Fur- 
ther, within linear response, the generalized susceptibil- 
ity G obtained from the slope of the EOS is related to the 
q = component of the equal-time correlation function, 
fc B TG(q = 0) = C(q = 0, t = 0) = / dr(Sp(r + x)<5p(x)) 
where 5p is the deviation from the mean density [T7]. We 
check whether analogous thermodynamic relations hold 
for the droplets taken as a subsystem in contact with the 
rest of the solid. We want to see how far these ther- 
modynamic considerations apply to our configurations of 
droplets, with the caveat that these averages are over 
a restricted ensemble and therefore does not affect the 
equilibrium thermodynamics of the bulk solid. 
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FIG. 5: (color-online) (a) Excess pressure, Ap c vs density p c of 
the clusters as a scatter plot at densities p — 0.88— 1.0 shown 
in different colors. Inspection of the individual droplets cor- 
responding to each of the colored dots shows that for Ap c > 
droplets tend to be string-like while Ap c < gives rise to com- 
pact droplets, (b) The radius of gyration R g of the droplets 
as a function of the number of particles in the droplets n c for 
several T — 0.35 (open squares), 0.40 (filled squares), 0.45 
(open circles) and 0.50 (filled circles). The data points for 
T > 0.35 have been each shifted by 1 to make them visible. 
Note that while R g ~ n c for Ap c < (i) it is linear for 
droplets with Ap c > (ii) (shown, in each case by red lines). 



data (Figjjjjb)) is consistent with Eq|3| however since 
the size of the clusters is not very large, it is difficult 
to probe the asymptotic behaviour. To obtain better 
statistics one needs to simulate yet larger systems for 
much longer times. 



B. Local thermodynamics of droplets 

The scatter diagram in Fig[5ja), suggests that there 
might be a thermodynamic interpretation of the local 



Pc 



FIG. 6: (color-online) (a) Scatter plot of p c vs p c of the 
droplets at T — 0.40 (red dots, the same data as in Figj5[a)) 
whose mean show shows a dependence similar to an "equa- 
tion of state" (green curve). The equilibrium equation of state 
of the solid is shown as a blue curve for comparison, (b) The 
scaled excess pressure p as a function of mean density p c for 
n c — 20 particle clusters for T — 0.35 — 1.0 in different colors. 
The data for T > 0.35 have each been shifted by 4 for an unit 
increase in T to make them visible. Note the prominent loop 
showing two distinct branches at low temperatures which dis- 
appears as T increases. We have obtained similar data for 
other values n c . 



For a fixed p and T for the solid, the density p c and 
pressure p c of the droplets is shown as a scatter plot in 
FigJ6^a). The averaged curve (green line in FigJ6]ja)) has 
a locus, p c vs p c which is distinct from the EOS of the 
equilibrium solid (blue line in FigJ6[a)). 

A plot of the scaled excess pressure p with the mean 
droplet density p c (Fig|6[b)) for fixed particle number 
n c shows a prominent non-linear feature akin to a van 
der Waals loop in the equilibrium pressure - density 
curve at a typical first order transition, say between gas 
and liquid. The two branches in the curves shown in 
Fig|6[b) correspond to compact (liquid-like) and string- 
like (glassy) droplets as discussed above. Integrating this 
pressure - density curve at fixed n c , gives the work done 
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FIG. 7: (color-online) Snapshots for T = 0.4 and p = 0.92 at 
0, 100, 150, 200, 250, 300, 350, 360, and 370 MD steps ((i)-(ix)) 
of a portion of our simulation cell showing the dynamics of 
two chosen droplets. The colors red and blue denotes +ve 
and — ve excess pressures respectively. Note that string-like 
droplets have Ap c > while relatively compact droplets have 
Ap c < 0. Note also the dissociation of a (red) string-like 
droplet into a (blue) compact and a smaller (red) string-like 
droplet ( (iv) — > (v)). 
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FIG. 8: The fluctuation -response ratio 7 shown as a function 
of n c for droplets with +ve (filled circles) and — ve excess pres- 
sure, p is independent of n c , with an intercept oc T. We verify 
this by plotting the intercept versus T for all the droplets (in- 
set) for T < 0.5. The filled and unfilled squares correspond 
to the + ve and — ve p branches respectively. 
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C. Density correlations 



by thermal fluctuations in creating a droplet of size n c . 
We find that at low temperatures, large excess pressure 
tends to convert compact droplets to stringy ones and 
vice versa (Figwl) over a characteristic relaxation time. 
Droplets are also seen to dissociate into distinct string- 
like and compact fragments with appropriate values for 
the excess pressure. Quite analogous to the familiar gas- 
liquid transition, this metastable "van der Waals loop " 
vanishes at higher temperatures beyond a "critical point" 
which exists somewhere in the range T = 0.9 — 1.0. 
Above this temperature the distinction between com- 
pact and string-like droplets disappears. We find that 
the generalized susceptibility G obtained from the slope 
of dp c /dp c is proportional to fluctuations of the density 
(Sp 2 ) = ((pc ~ Pc) 2 ) where (• • • ) denotes a time average. 
The fluctuation response ratio 7 = (5p 2 )n c / (p 2 c dp c /dp c ) 
should be independent of ji c , with an intercept which 
should be proportional to the temperature T. This is 
shown in FigjSlwhere 7 for the +ve and — ve p branches in 
Figpjjb) are shown separately using filled and open circles 
respectively. While 7 is indeed a constant, the intercept 
plotted versus T for all the droplets for the two branches 
are approximately linear in T. This is the metastable 
analogue of the equilibrium fluctuation-response relation 
discussed above, suggesting that the droplets are fluctua- 
tions from a metastable state describable by a free energy 
functional. 



We now study density correlations within each droplet 
to further characterize compact and string-like droplets. 
Note that these droplets have a finite lifetime t(ii c ), a 
plot of t vs. n c for fixed p and T is shown in Figj9j To 
obtain r, we have collected the times each constituent 
particle continues to belong to a droplet. As the droplet 
fluctuates, particles from the periphery continuously at- 
tach to and detach from the droplet leaving a set of par- 
ticles at the core intact. This is clear from the snapshots 
shown in Fig|4] We define r as the persistence time of 
these core particles. Our analysis for different p shows 
that the dense stringy clusters live longer. To obtain 
good statistics for the equal and unequal time density 
correlators, we therefore need to look at large and long- 



lived droplets. Fig 10 shows the equal time density cor- 
relations, g(r) . First of all g(r) in both branches show 
features associated with an amorphous structure with a 
first peak value which is much reduced from that of the 
full solid. Note that while the low density droplets are 
more liquid-like with smoothened peaks, the high den- 
sity droplets are glassy showing a prominent split sec- 
ond peak. We have also computed the self part of the 
van Hove correlation functions G s (t) = (p(0,0)p(0,t)). 
This is shown in Fig 11 While G s (t) relaxes exponen- 



tially for the low density, liquid-like droplets, the high 
density droplets show non-Debye relaxation. We show 
that as the lifetime of the droplets increase, a prominent 
/3-relaxation type plateau begins to develop. This, how- 
ever, gets cut-off by the finite lifetime of the droplets in 
the solid. 
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FIG. 9: Lifetime r vs n c for droplets at p = 0.92 and T = 0.4. 
Larger droplets survive longer. The error bars are obtained 
from the width of the distribution of r. The dashed straight 
line is a guide to the eye. 
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FIG. 10: Typical pair distribution function g(r) for particles 
in non-afnne droplets at p — 0.88 and T = 0.40, drawn sepa- 
rately for Ap c > (solid line) and Ap c < (dashed line) . The 
g(r) for the high density droplet has been shifted by 0.5 to 
make it visible. Both the g(r)'s show less crystalline structure 
than the solid (inset). 



IV. DISCUSSION AND CONCLUSIONS 

We have shown in this paper that the excitations of 
an equilibrium solid at high temperatures can be inter- 
preted as arising from a distribution of non-affine droplets 
whose mean size and lifetime increases as one approaches 
the liquid-solid phase boundary. These droplets are char- 
acterized by a density, internal pressure and shape. The 
shape of the droplets crossover from being compact to 
string-like as the density increases. The observed rela- 
tionship between the local pressure and density of the 
droplets in the form of an "equation of state" and the 
fluctuation-response relation of the local density strongly 
suggests that these non-afnne droplets arise as fluctua- 
tions from a metastable liquid or glass. Consistent with 
this we find that the high density droplets have a g(r) and 
van Hove function resembling that of glasses. Our main 
concern in this paper has been to provide evidence that 



10 u 



«J 10" 



o 



10' 







\. 03 10_1 
^"^"^^^ 10" ? 


V 


- \ ^ V * S *N. 


\^ t x 10 3 



20 



40 



60 



80 



100 



120 



t x 10 d 



FIG. 11: (a) The self part of the van-Hove function G s (t) as 
a function of time for the same droplets whose g(r) is given 
in Fig. |10| Note the prominent plateau for the high density 
droplets which is cut off by the finite lifetime of the droplet, 
(b) Development of the plateau in G s (t) with the lifetime of 
the droplets over which G s (t) is averaged (lines with mean r 
from 50 (bottom) to 200 (top) xl0~ 3 LJ units). 



droplet fluctuations of metastable liquid or glass exist in 
a crystalline solid and to characterize their shape and 
local thermodynamic parameters. How do these fluctua- 
tions influence the properties of solids? We hope to sys- 
tematically study and answer this question in the future. 
Some of the specific areas where the impact of droplet 
fluctuations may be observable are discussed below. 

For example, we wonder whether our results hint at 
the presence of a metastable liquid-glass critical point. 
While FigMb) is certainly suggestive, we must remember 
that the size of the non-affine droplets are typically small 
with n c ~ 100 even for the largest droplets. A careful 
finite size scaling analysis needs to be carried out in or- 
der to determine whether this feature survives for larger 
droplet sizes. The nature of the melting transition in the 
two-dimensional LJ system remains unclear in spite of 
being the subject of many investigations [5J fT5rf20"] over 
several decades. Early simulations obtained a first order 
liquid-solid and gas-solid transitions with prominent co- 
existence regions [51 [T5]. While at low temperatures near 
the triple point (T = 0.4), one obtains a first order melt- 
ing transitionfTH]. at higher temperatures (T ~ 1) one 
obtains a much reduced coexistence region jT9j with some 
characteristics of continuous [2 1J melting. At still higher 
temperatures, the melting transition appears to be un- 
equivocally driven by dislocation unbinding[20]. Exactly 
how and at what temperature one obtains this change in 
the nature of melting is, as yet, unknown. While we do 
not address this question in the present paper, we specu- 
late that droplet fluctuations may have a strong influence 
on the dynamics of melting. The presence of metastable 
critical points is known to crucially influence the dynam- 
ics of first order transitions, for eg. the important prob- 



lem of protein crvstallization [2"2l 123] . 

We expect these non-affine droplets to also play an 
important role in the rheology of solids under applied 
stresses; our preliminary work in this direction is consis- 
tent with this expectation. We would also like to enquire 
whether it is possible to observe these droplet excitations 
in a real experimental situation in two and three dimen- 
sions. Such excitations, if they exist, may be difficult to 
disentangle from the contributions coming from a den- 
sity of dislocations and grain boundaries. Perhaps direct 
visualization of droplet fluctuations in colloidal crystals 
is the best way to study these effects 24 . 
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